home *** CD-ROM | disk | FTP | other *** search
/ Software Vault: The Diamond Collection / The Diamond Collection (Software Vault)(Digital Impact).ISO / cdr44 / newmat08.zip / NL_EX.CPP < prev    next >
C/C++ Source or Header  |  1995-01-11  |  3KB  |  81 lines
  1. // This is an example of a non-linear least squares fit. The example
  2. // is from "Nonlinear estimation" by Gavin Ross (Springer,1990), p 63.
  3. // There are better ways of doing the fit in this case so this
  4. // example is just an example.
  5.  
  6. // The model is E(y) = a + b exp(-kx) and there are 6 data points.
  7.  
  8. #define WANT_STREAM
  9. #define WANT_MATH
  10. #include "newmatnl.h"
  11. #include "newmatio.h"
  12.  
  13. // first define the class describing the predictor function
  14.  
  15. class Model_3pe : public R1_Col_I_D
  16. {
  17.    ColumnVector x_values;         // the values of "x"
  18.    RowVector deriv;               // values of derivatives
  19. public:
  20.    Model_3pe(const ColumnVector& X_Values)
  21.       : x_values(X_Values) { deriv.ReDimension(3); }
  22.                                              // load X data
  23.    Real operator()(int);
  24.    Boolean IsValid() { return para(3)>0; }
  25.                                   // require "k" > 0
  26.    ReturnMatrix Derivatives() { return deriv; }
  27. };
  28.  
  29. Real Model_3pe::operator()(int i)
  30. {
  31.    Real a = para(1); Real b = para(2); Real k = para(3);
  32.    Real xvi = x_values(i);
  33.    Real e = exp(-k * xvi);
  34.    deriv(1) = 1.0;                    // calculate derivatives
  35.    deriv(2) = e;
  36.    deriv(3) = - b * e * xvi;
  37.    return a + b * e;                  // function value
  38. }
  39.  
  40. main()
  41. {
  42.    // Get the data
  43.    ColumnVector X(6);
  44.    ColumnVector Y(6);
  45.    X << 1   << 2   <<  3   <<  4   <<  6   <<  8;
  46.    Y << 3.2 << 7.9 << 11.1 << 14.5 << 16.7 << 18.3;
  47.  
  48.  
  49.    // Do the fit
  50.    Model_3pe model(X);                // the model object
  51.    NonLinearLeastSquares NLLS(model); // the non-linear least squares object
  52.    ColumnVector Para(3);              // for the parameters
  53.    Para << 9 << -6 << .5;             // trial values of parameters
  54.    cout << "Fitting parameters\n";
  55.    NLLS.Fit(Y,Para);                  // do the fit
  56.  
  57.    // Inspect the results
  58.    ColumnVector SE;                   // for the standard errors
  59.    NLLS.GetStandardErrors(SE);
  60.    cout << "\n\nEstimates and standard errors\n" <<
  61.       setw(10) << setprecision(2) << (Para | SE) << endl;
  62.    Real ResidualSD = sqrt(NLLS.ResidualVariance());
  63.    cout << "\nResidual s.d. = " << setw(10) << setprecision(2) <<
  64.       ResidualSD << endl;
  65.    SymmetricMatrix Correlations;
  66.    NLLS.GetCorrelations(Correlations);
  67.    cout << "\nCorrelationMatrix\n" <<
  68.       setw(10) << setprecision(2) << Correlations << endl;
  69.    ColumnVector Residuals;
  70.    NLLS.GetResiduals(Residuals);
  71.    DiagonalMatrix Hat;
  72.    NLLS.GetHatDiagonal(Hat);
  73.    cout << "\nX, Y, Residual, Hat\n" << setw(10) << setprecision(2) <<
  74.       (X | Y | Residuals | Hat.AsColumn()) << endl;
  75.    // recover var/cov matrix
  76.    SymmetricMatrix D;
  77.    D << SE.AsDiagonal() * Correlations * SE.AsDiagonal();
  78.    cout << "\nVar/cov\n" << setw(14) << setprecision(4) << D << endl;
  79.    return 0;
  80. }
  81.